### Syntax Democratic Deconsolidation
### APC-Analysis on democratic support in Europe


#Get all packages with the version used for this analysis
renv::restore()

#Load all packages required for analysis

if (!require("pacman")) install.packages("pacman")
library(pacman)

p_load(tidyverse, haven,  pastecs, mgcv, itsadug, lme4, broom, sjPlot, gridExtra, grid, lattice, merTools, rlist, texreg, ggrepel, ragg)

#Set working directory
wd <- here::here()
setwd(wd)

###Take care that the subfolder data exists and the datasets "dataEVS_WVS.dta" and "ZA7500_v2-0-0.dta" exist!


#Create necessary sub-directories for plots, appendix and Shiny App
dir.create("plots",showWarnings = F)
dir.create("ShinyPlots", showWarnings = F)
dir.create("ShinyData", showWarnings = F)

wdPlots <- paste0(wd, "/plots/")
wdShinyPlots <- paste0(wd, "/ShinyPlots/")
wdShinyData <- paste0(wd, "/ShinyData/")


##Load data
wvs_evs <- read_dta("data/dataEVS_WVS.dta")
w5 <- read_dta("data/ZA7500_v2-0-0.dta")


#Create notin function
"%notin%" <- Negate("%in%")



##### Country selection

##List of all countries that will be reported in the appendix (maintext + appendixonly)
country_allcountries <-  c("Albania", "Austria", "Bulgaria", "Croatia", "Czech Republic", "Denmark", "Estonia", "Finland", "France", "Germany", "Hungary", "Iceland", "Italy", "Lithuania", "Netherlands", "Norway", "Poland", "Romania", "Slovakia",       "Slovenia", "Spain", "Sweden", "Switzerland", "United Kingdom")

##List of countries that will be reported !only! in the appendix   
country_appendixonly  <- c("Albania", "Bulgaria", "Croatia", "Czech Republic", "Estonia", "Romania")

#List of countries, where all indicators are available
countrySub <- c("Bulgaria", "Estonia", "Finland", "France", "Germany", "Hungary", "Italy", "Netherlands", "Norway", "Poland", "Romania", "Slovenia", "Spain", "Sweden", "Switzerland", "United Kingdom")

##List of countries that will be reported in the main text
country_maintext <- setdiff(country_allcountries, country_appendixonly)



#####Data Transformation#####


##Recoding of pre-2017 survey data
wvs_evs_all <- wvs_evs %>% 
  mutate(country = as.factor(S003),
         country = fct_recode(country, "Albania" = "8", "Azerbaijan" = "31", "Austria" = "40", "Armenia" = "51", "Bulgaria" = "100", "Belarus" = "112", "Croatia" = "191", "Czech Republic" = "203", "Denmark" = "208", "Estonia" = "233", "Finland" = "246", "France" = "250", "Georgia" = "268", "Germany" = "276", "Hungary" = "348", "Iceland" = "352", "Italy" = "380", "Lithuania" = "440", "Netherlands" = "528", "Norway" = "578", "Poland" = "616", "Romania" = "642", "Russian Federation" = "643", "Serbia" = "688", "Slovakia" = "703", "Slovenia" = "705", "Spain" = "724", "Sweden" = "752", "Switzerland" = "756", "United Kingdom" = "826"),
         wavesWVS = as.factor(S002),
         wavesWVS = fct_recode(wavesWVS, "WVS2" = "2", "WVS3" = "3", "WVS4" = "4", "WVS5" = "5", "WVS6" = "6"),
         wavesEVS = as.factor(S002EVS),
         wavesEVS = fct_recode(wavesEVS, "EVS1" = "1", "EVS2" = "2", "EVS3" = "3", "EVS4" = "4", NULL = "-4")) %>% 
  mutate(waves = ifelse(!is.na(wavesWVS), levels(wavesWVS)[wavesWVS], levels(wavesEVS)[wavesEVS]),
         waves = as.factor(waves),
         weights = as.numeric(S017)) %>% 
  mutate(polImportant = if_else(as.numeric(A004) < 1, NaN, as.numeric(A004)),
         polImportant = (polImportant - 4) / (-3), #Standardized
         authoritarism = if_else(as.numeric(E018) < 1, NaN, as.numeric(E018)),
         authoritarism = (authoritarism - 2) * (-1), #Good thing at 1, bad thing at -1
         polint = if_else(as.numeric(E023) < 1, NaN, as.numeric(E023)),
         polint = (polint - 4) / (-3), #Standardized, do not use EVS 1981 (not included in our sample, so no problem)
         system_democracy = if_else(as.numeric(E117) < 1, NaN, as.numeric(E117)),
         system_democracy = (system_democracy - 4) / (-3),#Standardized
         system_leader = if_else(as.numeric(E114) < 1, NaN, as.numeric(E114)),
         system_leader = (system_leader - 4) / (-3),#Standardized
         system_experts = if_else(as.numeric(E115) < 1, NaN, as.numeric(E115)),
         system_experts = (system_experts - 4) / (-3),#Standardized
         system_army = if_else(as.numeric(E116) < 1, NaN, as.numeric(E116)),
         system_army = (system_army - 4) / (-3),#Standardized
         democracy_freeElection = if_else(as.numeric(E226) < 1, NaN, as.numeric(E226)),
         democracy_freeElection = (democracy_freeElection - 1) / (9),#Standardized
         democracy_ArmyTakeOver = if_else(as.numeric(E228) < 1, NaN, as.numeric(E228)),
         democracy_ArmyTakeOver = (democracy_ArmyTakeOver - 1) / (9),#Standardized
         democracy_protectLiberty = if_else(as.numeric(E229) < 1, NaN, as.numeric(E229)),
         democracy_protectLiberty = (democracy_protectLiberty - 1) / (9),#Standardized
         DemImportant = if_else(as.numeric(E235) < 1, NaN, as.numeric(E235)),
         DemImportant = (DemImportant - 1) / (9),#Standardized
         trustParliament = if_else(as.numeric(E069_07) < 1, NaN, as.numeric(E069_07)),
         trustParliament = (trustParliament - 4) / (-3),#Standardized
         trustJustice = if_else(as.numeric(E069_17) < 1, NaN, as.numeric(E069_17)),
         trustJustice = (trustJustice - 4) / (-3),#Standardized
         trustCivilService = if_else(as.numeric(E069_08) < 1, NaN, as.numeric(E069_08)),
         trustCivilService = (trustCivilService - 4) / (-3),#Standardized
         
         trustInstSmall = (trustParliament + trustJustice + trustCivilService) / 3, #Additive Index
         
         sex = if_else(as.numeric(X001) < 1, NaN, as.numeric(X001)),
         sex = sex - 1,
         birthyear = if_else(as.numeric(X002) < 1880, NaN, as.numeric(X002)), #Remove extremely unreasonable respondent ages (older than 110?)
         period = as.numeric(S020)
  ) %>% 
  filter(country %in% c("Austria", "Germany", "Iceland", "Netherlands", "Spain", "Switzerland", "Slovenia", "Bulgaria", "Croatia", "Czech Republic", "Poland", "Slovakia", "Denmark", "Finland", "France", "United Kingdom", "Italy", "Sweden", "Hungary", "Lithuania", "Norway", "Estonia", "Romania", "Albania")) %>% 
  dplyr::select(country, waves, weights, polImportant, authoritarism, polint, system_democracy, system_leader, system_experts, system_army, democracy_freeElection, democracy_ArmyTakeOver, democracy_protectLiberty, DemImportant, trustParliament, trustJustice, trustCivilService, trustInstSmall, sex, birthyear, period)



##Recoding of 2017 survey data of the EVS
w5_all <- w5 %>% 
  mutate(country = as.factor(country),
         country = fct_recode(country, "Albania" = "8", "Azerbaijan" = "31", "Austria" = "40", "Armenia" = "51", "Bulgaria" = "100", "Belarus" = "112", "Croatia" = "191", "Czech Republic" = "203", "Denmark" = "208", "Estonia" = "233", "Finland" = "246", "France" = "250", "Georgia" = "268", "Germany" = "276", "Hungary" = "348", "Iceland" = "352", "Italy" = "380", "Lithuania" = "440", "Netherlands" = "528", "Norway" = "578", "Poland" = "616", "Romania" = "642", "Russian Federation" = "643", "Serbia" = "688", "Slovakia" = "703", "Slovenia" = "705", "Spain" = "724", "Sweden" = "752", "Switzerland" = "756", "United Kingdom" = "826"),
         wavesEVS = "EVS5") %>% 
  filter(country %in% c("Austria", "Germany", "Iceland", "Netherlands", "Spain", "Switzerland", "Slovenia", "Bulgaria", "Croatia", "Czech Republic", "Poland", "Slovakia", "Denmark", "Finland", "France", "United Kingdom", "Italy", "Sweden", "Hungary", "Lithuania", "Norway", "Estonia", "Romania", "Albania")) %>% 
  mutate(waves = wavesEVS,
         waves = as.factor(waves),
         weights = gweight) %>% 
  mutate(polImportant = if_else(as.numeric(v5) < 1, NaN, as.numeric(v5)),
         polImportant = (polImportant - 4) / (-3), #Standardized
         authoritarism = if_else(as.numeric(v114) < 1, NaN, as.numeric(v114)),
         authoritarism = (authoritarism - 2) * (-1), #Good thing at 1, bad thing at -1
         polint = if_else(as.numeric(v97) < 1, NaN, as.numeric(v97)),
         polint = (polint - 4) / (-3), #Standardized, do not use EVS 1981 (not included in our sample, so no problem)
         system_democracy = if_else(as.numeric(v148) < 1, NaN, as.numeric(v148)),
         system_democracy = (system_democracy - 4) / (-3),#Standardized
         system_leader = if_else(as.numeric(v145) < 1, NaN, as.numeric(v145)),
         system_leader = (system_leader - 4) / (-3),#Standardized
         system_experts = if_else(as.numeric(v146) < 1, NaN, as.numeric(v146)),
         system_experts = (system_experts - 4) / (-3),#Standardized
         system_army = if_else(as.numeric(v147) < 1, NaN, as.numeric(v147)),
         system_army = (system_army - 4) / (-3),#Standardized
         democracy_freeElection = if_else(as.numeric(v135) < 1, NaN, as.numeric(v135)),
         democracy_freeElection = (democracy_freeElection - 1) / (9),#Standardized
         democracy_ArmyTakeOver = if_else(as.numeric(v137) < 1, NaN, as.numeric(v137)),
         democracy_ArmyTakeOver = (democracy_ArmyTakeOver - 1) / (9),#Standardized
         democracy_protectLiberty = if_else(as.numeric(v138) < 1, NaN, as.numeric(v138)),
         democracy_protectLiberty = (democracy_protectLiberty - 1) / (9),#Standardized
         DemImportant = if_else(as.numeric(v142) < 1, NaN, as.numeric(v142)),
         DemImportant = (DemImportant - 1) / (9),#Standardized
         trustParliament = if_else(as.numeric(v121) < 1, NaN, as.numeric(v121)),
         trustParliament = (trustParliament - 4) / (-3),#Standardized
         trustJustice = if_else(as.numeric(v127) < 1, NaN, as.numeric(v127)),
         trustJustice = (trustJustice - 4) / (-3),#Standardized
         trustCivilService = if_else(as.numeric(v122) < 1, NaN, as.numeric(v122)),
         trustCivilService = (trustCivilService - 4) / (-3),#Standardized
         
         trustInstSmall = (trustParliament + trustJustice + trustCivilService) / 3, #Additive Index
         
         sex = if_else(as.numeric(v225) < 1, NaN, as.numeric(v225)),
         sex = sex - 1,
         birthyear = if_else(as.numeric(v226) < 1880, NaN, as.numeric(v226)), #Remove extremely unreasonable respondent ages (older than 110?)
         period = 2017,
         period = as.numeric(period)
  ) %>% 
  mutate(weights = pmin(pmax(weights, quantile(weights, .025)), 
                        quantile(weights, .975))) %>%
  dplyr::select(country, waves, weights, polImportant, authoritarism, polint, system_democracy, system_leader, system_experts, system_army, democracy_freeElection, democracy_ArmyTakeOver, democracy_protectLiberty, DemImportant, trustParliament, trustJustice, trustCivilService, trustInstSmall, sex, birthyear, period)



#Bind data
datComplete <- wvs_evs_all %>% 
  bind_rows(w5_all) %>% 
  mutate(
    age = period - birthyear
  )

#Create age and cohort groups
datComplete <- datComplete %>% 
  mutate(born_adult = birthyear + 18) %>% 
  mutate(age_grp = ifelse(age >= 15 & age <20, 1, ifelse(
    age >= 20 & age < 30,  2, ifelse(
      age >= 30 & age < 40, 3, ifelse(
        age >= 40 & age < 50,  4, ifelse(
          age >= 50 & age < 60,  5, ifelse(
            age >= 60 & age < 70,  6, ifelse(
              age >= 70 & age < 80,  7,  8
            )
          )
        )
      )
    )
  )), 
  age_grp_centered = age_grp - mean(age_grp, na.rm = T),
  born_adult_centered = born_adult - mean(born_adult, na.rm = T),
  age_grp1 = ifelse(age >= 15 & age < 30, 1,0),     ##Age groups for GAM model
  age_grp2 = ifelse(age >= 30 & age < 60, 1, 0),
  period.f = as.factor(period),
  waves.f = as.factor(waves)
  )
#Remove unreasonably young and old respondents

datComplete <- datComplete %>% 
  filter(age >= 16 & age <=95) %>% 
  mutate(age_cohort_five = findInterval(born_adult, c(seq(1925, 2010, by = 5))))  ##cohort groups for rHAPC


###Select the correct waves and countries for each analysis
datCompleteNew <- datComplete %>% 
  gather(dv_name, dv_value, c(polImportant:trustInstSmall)) %>% 
  filter(waves %in% c("EVS1","EVS2","EVS3","EVS4","EVS5", "WVS5", "WVS6")) %>% 
  filter(waves %in% c("EVS1","EVS2","EVS3","EVS4","EVS5") | 
           ((waves %in% c("WVS5")) & country %in% countrySub) | 
           ((waves %in% c("WVS6")) & country %in% countrySub)) %>% 
  dplyr::select(country, waves, sex, birthyear, period, age, born_adult, age_grp, age_grp_centered, born_adult_centered, age_grp1, age_grp2, period.f, waves.f, age_cohort_five, dv_name, dv_value) %>% 
  drop_na()


#####Data Transformation - only Interviewed considered in Wave 5#####


#####Recoding of pre-2017 survey data##########
####Directly downloaded from the EVS homepage
wvs_evs_interviewed <- wvs_evs %>% 
  mutate(country = as.factor(S003),
         country = fct_recode(country, "Albania" = "8", "Azerbaijan" = "31", "Austria" = "40", "Armenia" = "51", "Bulgaria" = "100", "Belarus" = "112", "Croatia" = "191", "Czech Republic" = "203", "Denmark" = "208", "Estonia" = "233", "Finland" = "246", "France" = "250", "Georgia" = "268", "Germany" = "276", "Hungary" = "348", "Iceland" = "352", "Italy" = "380", "Lithuania" = "440", "Netherlands" = "528", "Norway" = "578", "Poland" = "616", "Romania" = "642", "Russian Federation" = "643", "Serbia" = "688", "Slovakia" = "703", "Slovenia" = "705", "Spain" = "724", "Sweden" = "752", "Switzerland" = "756", "United Kingdom" = "826"),
         wavesWVS = as.factor(S002),
         wavesWVS = fct_recode(wavesWVS, "WVS2" = "2", "WVS3" = "3", "WVS4" = "4", "WVS5" = "5", "WVS6" = "6"),
         wavesEVS = as.factor(S002EVS),
         wavesEVS = fct_recode(wavesEVS, "EVS1" = "1", "EVS2" = "2", "EVS3" = "3", "EVS4" = "4", NULL = "-4")) %>% 
  mutate(waves = ifelse(!is.na(wavesWVS), levels(wavesWVS)[wavesWVS], levels(wavesEVS)[wavesEVS]),
         waves = as.factor(waves),
         weights = as.numeric(S017)) %>% 
  mutate(polImportant = if_else(as.numeric(A004) < 1, NaN, as.numeric(A004)),
         polImportant = (polImportant - 4) / (-3), #Standardized
         authoritarism = if_else(as.numeric(E018) < 1, NaN, as.numeric(E018)),
         authoritarism = (authoritarism - 2) * (-1), #Good thing at 1, bad thing at -1
         polint = if_else(as.numeric(E023) < 1, NaN, as.numeric(E023)),
         polint = (polint - 4) / (-3), #Standardized, do not use EVS 1981 (not included in our sample, so no problem)
         system_democracy = if_else(as.numeric(E117) < 1, NaN, as.numeric(E117)),
         system_democracy = (system_democracy - 4) / (-3),#Standardized
         system_leader = if_else(as.numeric(E114) < 1, NaN, as.numeric(E114)),
         system_leader = (system_leader - 4) / (-3),#Standardized
         system_experts = if_else(as.numeric(E115) < 1, NaN, as.numeric(E115)),
         system_experts = (system_experts - 4) / (-3),#Standardized
         system_army = if_else(as.numeric(E116) < 1, NaN, as.numeric(E116)),
         system_army = (system_army - 4) / (-3),#Standardized
         democracy_freeElection = if_else(as.numeric(E226) < 1, NaN, as.numeric(E226)),
         democracy_freeElection = (democracy_freeElection - 1) / (9),#Standardized
         democracy_ArmyTakeOver = if_else(as.numeric(E228) < 1, NaN, as.numeric(E228)),
         democracy_ArmyTakeOver = (democracy_ArmyTakeOver - 1) / (9),#Standardized
         democracy_protectLiberty = if_else(as.numeric(E229) < 1, NaN, as.numeric(E229)),
         democracy_protectLiberty = (democracy_protectLiberty - 1) / (9),#Standardized
         DemImportant = if_else(as.numeric(E235) < 1, NaN, as.numeric(E235)),
         DemImportant = (DemImportant - 1) / (9),#Standardized
         trustParliament = if_else(as.numeric(E069_07) < 1, NaN, as.numeric(E069_07)),
         trustParliament = (trustParliament - 4) / (-3),#Standardized
         trustJustice = if_else(as.numeric(E069_17) < 1, NaN, as.numeric(E069_17)),
         trustJustice = (trustJustice - 4) / (-3),#Standardized
         trustCivilService = if_else(as.numeric(E069_08) < 1, NaN, as.numeric(E069_08)),
         trustCivilService = (trustCivilService - 4) / (-3),#Standardized
         
         trustInstSmall = (trustParliament + trustJustice + trustCivilService) / 3, #Additive Index
         
         sex = if_else(as.numeric(X001) < 1, NaN, as.numeric(X001)),
         sex = sex - 1,
         birthyear = if_else(as.numeric(X002) < 1880, NaN, as.numeric(X002)), #Remove extremely unreasonable respondent ages (older than 110?)
         period = as.numeric(S020)
  ) %>% 
  filter(country %in% c("Austria", "Germany", "Iceland", "Netherlands", "Spain", "Switzerland", "Slovenia", "Bulgaria", "Croatia", "Czech Republic", "Poland", "Slovakia", "Denmark", "Finland", "France", "United Kingdom", "Italy", "Sweden", "Hungary", "Lithuania", "Norway", "Estonia", "Romania", "Albania")) %>% 
  dplyr::select(country, waves, weights, polImportant, authoritarism, polint, system_democracy, system_leader, system_experts, system_army, democracy_freeElection, democracy_ArmyTakeOver, democracy_protectLiberty, DemImportant, trustParliament, trustJustice, trustCivilService, trustInstSmall, sex, birthyear, period)



#######Recoding of 2017 survey data of the EVS##########################
w5_interviewed <- w5 %>% 
  mutate(country = as.factor(country),
         country = fct_recode(country, "Albania" = "8", "Azerbaijan" = "31", "Austria" = "40", "Armenia" = "51", "Bulgaria" = "100", "Belarus" = "112", "Croatia" = "191", "Czech Republic" = "203", "Denmark" = "208", "Estonia" = "233", "Finland" = "246", "France" = "250", "Georgia" = "268", "Germany" = "276", "Hungary" = "348", "Iceland" = "352", "Italy" = "380", "Lithuania" = "440", "Netherlands" = "528", "Norway" = "578", "Poland" = "616", "Romania" = "642", "Russian Federation" = "643", "Serbia" = "688", "Slovakia" = "703", "Slovenia" = "705", "Spain" = "724", "Sweden" = "752", "Switzerland" = "756", "United Kingdom" = "826"),
         wavesEVS = "EVS5") %>% 
  filter(country %in% c("Austria", "Germany", "Iceland", "Netherlands", "Spain", "Switzerland", "Slovenia", "Bulgaria", "Croatia", "Czech Republic", "Poland", "Slovakia", "Denmark", "Finland", "France", "United Kingdom", "Italy", "Sweden", "Hungary", "Lithuania", "Norway", "Estonia", "Romania", "Albania")) %>% 
  mutate(waves = wavesEVS,
         waves = as.factor(waves),
         weights = gweight) %>% 
  mutate(polImportant = if_else(as.numeric(v5) < 1, NaN, as.numeric(v5)),
         polImportant = (polImportant - 4) / (-3), #Standardized
         authoritarism = if_else(as.numeric(v114) < 1, NaN, as.numeric(v114)),
         authoritarism = (authoritarism - 2) * (-1), #Good thing at 1, bad thing at -1
         polint = if_else(as.numeric(v97) < 1, NaN, as.numeric(v97)),
         polint = (polint - 4) / (-3), #Standardized, do not use EVS 1981 (not included in our sample, so no problem)
         system_democracy = if_else(as.numeric(v148) < 1, NaN, as.numeric(v148)),
         system_democracy = (system_democracy - 4) / (-3),#Standardized
         system_leader = if_else(as.numeric(v145) < 1, NaN, as.numeric(v145)),
         system_leader = (system_leader - 4) / (-3),#Standardized
         system_experts = if_else(as.numeric(v146) < 1, NaN, as.numeric(v146)),
         system_experts = (system_experts - 4) / (-3),#Standardized
         system_army = if_else(as.numeric(v147) < 1, NaN, as.numeric(v147)),
         system_army = (system_army - 4) / (-3),#Standardized
         democracy_freeElection = if_else(as.numeric(v135) < 1, NaN, as.numeric(v135)),
         democracy_freeElection = (democracy_freeElection - 1) / (9),#Standardized
         democracy_ArmyTakeOver = if_else(as.numeric(v137) < 1, NaN, as.numeric(v137)),
         democracy_ArmyTakeOver = (democracy_ArmyTakeOver - 1) / (9),#Standardized
         democracy_protectLiberty = if_else(as.numeric(v138) < 1, NaN, as.numeric(v138)),
         democracy_protectLiberty = (democracy_protectLiberty - 1) / (9),#Standardized
         DemImportant = if_else(as.numeric(v142) < 1, NaN, as.numeric(v142)),
         DemImportant = (DemImportant - 1) / (9),#Standardized
         trustParliament = if_else(as.numeric(v121) < 1, NaN, as.numeric(v121)),
         trustParliament = (trustParliament - 4) / (-3),#Standardized
         trustJustice = if_else(as.numeric(v127) < 1, NaN, as.numeric(v127)),
         trustJustice = (trustJustice - 4) / (-3),#Standardized
         trustCivilService = if_else(as.numeric(v122) < 1, NaN, as.numeric(v122)),
         trustCivilService = (trustCivilService - 4) / (-3),#Standardized
         
         trustInstSmall = (trustParliament + trustJustice + trustCivilService) / 3, #Additive Index
         
         sex = if_else(as.numeric(v225) < 1, NaN, as.numeric(v225)),
         sex = sex - 1,
         birthyear = if_else(as.numeric(v226) < 1880, NaN, as.numeric(v226)), #Remove extremely unreasonable respondent ages (older than 110?)
         period = 2017,
         period = as.numeric(period)
  ) %>% 
  mutate(weights = pmin(pmax(weights, quantile(weights, .025)), 
                        quantile(weights, .975))) %>%
  filter(mm_select_sample == 1) %>% 
  dplyr::select(country, waves, weights, polImportant, authoritarism, polint, system_democracy, system_leader, system_experts, system_army, democracy_freeElection, democracy_ArmyTakeOver, democracy_protectLiberty, DemImportant, trustParliament, trustJustice, trustCivilService, trustInstSmall, sex, birthyear, period)




#Bind data
datComplete_interviewed <- wvs_evs_interviewed %>% 
  bind_rows(w5_interviewed) %>% 
  mutate(
    age = period - birthyear
  )

#Create age and cohort groups
datComplete_interviewed <- datComplete_interviewed %>% 
  mutate(born_adult = birthyear + 18) %>% 
  mutate(age_grp = ifelse(age >= 15 & age <20, 1, ifelse(
    age >= 20 & age < 30,  2, ifelse(
      age >= 30 & age < 40, 3, ifelse(
        age >= 40 & age < 50,  4, ifelse(
          age >= 50 & age < 60,  5, ifelse(
            age >= 60 & age < 70,  6, ifelse(
              age >= 70 & age < 80,  7,  8
            )
          )
        )
      )
    )
  )), 
  age_grp_centered = age_grp - mean(age_grp, na.rm = T),
  born_adult_centered = born_adult - mean(born_adult, na.rm = T),
  age_grp1 = ifelse(age >= 15 & age < 30, 1,0),     ##Age groups for GAM model
  age_grp2 = ifelse(age >= 30 & age < 60, 1, 0),
  period.f = as.factor(period),
  waves.f = as.factor(waves)
  )
#Remove unreasonably young and old respondents

datComplete_interviewed <- datComplete_interviewed %>% 
  filter(age >= 16 & age <=95) %>% 
  mutate(age_cohort_five = findInterval(born_adult, c(seq(1925, 2010, by = 5))))  ##cohort groups for rHAPC


###Select the correct waves and countries for each analysis
datCompleteNew_interviewed <- datComplete_interviewed %>% 
  gather(dv_name, dv_value, c(polImportant:trustInstSmall)) %>% 
  filter(waves %in% c("EVS1","EVS2","EVS3","EVS4","EVS5", "WVS5", "WVS6")) %>% 
  filter(waves %in% c("EVS1","EVS2","EVS3","EVS4","EVS5") | 
           ((waves %in% c("WVS5")) & country %in% countrySub) | 
           ((waves %in% c("WVS6")) & country %in% countrySub)) %>% 
  dplyr::select(country, waves, sex, birthyear, period, age, born_adult, age_grp, age_grp_centered, born_adult_centered, age_grp1, age_grp2, period.f, waves.f, age_cohort_five, dv_name, dv_value) %>% 
  drop_na()


#####Mean Trends over Time: Broad Country Sample for Main text #####
#Figure 1
descPlots_maintext <- datCompleteNew_interviewed %>% 
  filter(country %in% country_maintext) %>% 
  filter(waves %in% c("EVS1","EVS2","EVS3","EVS4","EVS5")) %>% 
  filter(!(dv_name == "democracy_freeElection" | dv_name == "democracy_ArmyTakeOver" | dv_name == "democracy_protectLiberty" | dv_name == "DemImportant")) %>% 
  group_by(dv_name, country, waves) %>% 
  summarize(dv_value = mean(dv_value)) %>% 
  filter(dv_name %in% c("system_democracy", "system_leader", "trustInstSmall")) %>% 
  ungroup() %>% 
  mutate(dv_name = factor(dv_name, levels=c('system_democracy', 'trustInstSmall', 'system_leader'))) %>% 
  mutate(dv_name = recode(dv_name, system_democracy = "Regime: Democracy", trustInstSmall = "Trust Institutions", system_leader = "Regime: Strong Leader")) %>% 
  mutate(waves = factor(waves)) %>% 
  mutate(waves = fct_relevel(waves, "EVS5", "EVS4", "EVS3", "EVS2", "EVS1")) %>% 
  mutate(waves = fct_recode(waves, "2018" = "EVS5", "2008" = "EVS4", "1999" = "EVS3", "1990" = "EVS2", "1981" = "EVS1")) %>%
  ggplot(aes(x = waves, y = dv_value)) +
  geom_segment( aes(x=waves, xend=waves, y=0, yend=dv_value)) +
  geom_point(size=2, alpha=0.6) +
  ylim(0, 1) +
  coord_flip() + 
  theme_light() +
  theme(text = element_text(size=8), axis.text = element_text(size = 6), strip.text.y = element_text(angle=0)) + 
  facet_grid(country~dv_name) + 
  labs(x = "", y = "")

tiff(file=paste0(wdPlots, "figure1_trend_onlyInterviewed.tiff"), width = 2800, height = 3500, res = 400)

print(descPlots_maintext)

dev.off()


#Get ggplot2 colors
gg_color_hue <- function(n) {
  hues = seq(15, 375, length = n + 1)
  hcl(h = hues, l = 65, c = 100)[1:n]
}

cols = gg_color_hue(2)

#Subsample by items

countrySub_DemImportant <- countrySub

countrySub_Essentials <- countrySub[!countrySub %in% "Italy"]

#Run Primary models
source("02_APC_Analysis_GAM.R")
source("03_APC_Analysis_onlyInterviewed_GAM.R")

#Run secondary models
source("04_APC_Analysis_HAPC.R")
source("05_APC_Analysis_onlyInterviewed_HAPC.R")
